---
title: "Multi-Level Conceptual Replication of Fulmer et al. (2010): Table"
author: "Matthew Samson"
date: "19 January 2015"
output: html_document
---

Install and require appropriate packages.

```{r}

install.packages(c("psych","nlme", "reshape", "ICC", "plyr", "data.table"))

require("psych"); require("nlme"); require("reshape"); require("ICC"); require("plyr"); require("data.table")

```

Set working directory and load data. 

```{r, echo=TRUE}

#Set Working Directory (insert appropriate directory below)
setwd("~/mjs268/2_Direct Replication of Fulmer et al. (2010)/4_Data for Analysis")

#Load Data
ics2001 <- read.csv("3_International College Survey 2001 Final.csv")

```

Adjust for reverse coding and aggregate variables in ics2001

```{r}

#Reverse coding
ics2001$ext2 <- 6 - ics2001$ext2
ics2001$ext3 <- 6 - ics2001$ext3
ics2001$ext6 <- 6 - ics2001$ext6

#Average Extroversion
ics2001$ext_tot <- (ics2001$ext1 + ics2001$ext2 + ics2001$ext3 + ics2001$ext4 + ics2001$ext5 + ics2001$ext6)/6
#Average Subjective Well-Being
ics2001$swl_tot <- (ics2001$swl1 + ics2001$swl2 + ics2001$swl3 + ics2001$swl4 + ics2001$swl5)/5

```

Create temp file for all L1 variables (which are from ics2001) used in analysis

```{r}

TEMP1 <- (ics2001[,c("country", "acq", "ext_tot", "swl_tot", "ext1", "ext2", "ext3", "ext4", "ext5", "ext6", "swl1", "swl2", "swl3", "swl4", "swl5")])

#Omit participants flagged by Diener et al. (2001) as acquiescent. 
TEMP1 <- subset(TEMP1, acq == "0")

```

Create group means and SD's for EXT and SWL

```{r}

#Means
country.means <- data.table(TEMP1)
country.means <- data.frame(country.means[,list(country.ext_tot=mean(ext_tot,na.rm=T),country.swl_tot=mean(swl_tot,na.rm=T)),by=country])

#Rename
country.means <- rename(country.means, c("country.ext_tot"="EXT_M", "country.swl_tot"="SWL_M"))

#SD's
country.sds <- data.table(TEMP1)
country.sds <- data.frame(country.sds[,list(country.ext_tot=sd(ext_tot,na.rm=T),country.swl_tot=sd(swl_tot,na.rm=T)),by=country])

#Rename
country.sds <- rename(country.sds, c("country.ext_tot"="EXT_SD", "country.swl_tot"="SWL_SD"))

#Merge
country.means_sds <- merge(country.means,country.sds,by="country")
TEMP1 <- merge(TEMP1,country.means_sds,by="country")

```

Delete construct scores that are +/- 3 SD from country mean for EXT and SWL.

```{r}

#Compute Z-Scores
TEMP1$z_ext <- (TEMP1$ext_tot-TEMP1$EXT_M)/TEMP1$EXT_SD
TEMP1$z_swl <- (TEMP1$swl_tot-TEMP1$SWL_M)/TEMP1$SWL_SD

#Delete outliers 
TEMP1$ext_tot[TEMP1$z_ext>3 & !is.na(TEMP1$ext_tot)] <- NA
TEMP1$ext_tot[TEMP1$z_ext<(-3) & !is.na(TEMP1$ext_tot)] <- NA
TEMP1$swl_tot[TEMP1$z_swl>3 & !is.na(TEMP1$swl_tot)] <- NA
TEMP1$swl_tot[TEMP1$z_swl<(-3) & !is.na(TEMP1$swl_tot)] <- NA

```

Create group means and SD's for EXT and SWL **without outliers**

```{r}

#Delete now irrelevant variables from data set
TEMP1 <- na.omit(TEMP1[,c("country", "ext_tot", "swl_tot", "ext1", "ext2", "ext3", "ext4", "ext5", "ext6", "swl1", "swl2", "swl3", "swl4", "swl5")])

#Means
country.means <- data.table(TEMP1)
country.means <- data.frame(country.means[,list(country.ext_tot=mean(ext_tot,na.rm=T),country.swl_tot=mean(swl_tot,na.rm=T)),by=country])

#Rename
country.means <- rename(country.means, c("country.ext_tot"="EXT_M", "country.swl_tot"="SWL_M"))

#SD's
country.sds <- data.table(TEMP1)
country.sds <- data.frame(country.sds[,list(country.ext_tot=sd(ext_tot,na.rm=T),country.swl_tot=sd(swl_tot,na.rm=T)),by=country])

#Rename
country.sds <- rename(country.sds, c("country.ext_tot"="EXT_SD", "country.swl_tot"="SWL_SD"))

#Merge
country.means_sds <- merge(country.means,country.sds,by="country")

```

Internal consistency of extraversion scale by country.

```{r}

country_vector <- c()
ext_alpha_vector <- c()

for (temp.country in levels(TEMP1$country)) {
  extalpha <- subset(TEMP1, country==temp.country, ext1:ext6)
  a<-alpha(extalpha, keys = c(1,1,1,1,1,1))
  alpha_output<-as.numeric(a$total[1])
  country_vector<-append(country_vector, temp.country)
  ext_alpha_vector<-append(ext_alpha_vector,alpha_output)
  }

extroversion_alphas <- cbind(country_vector, ext_alpha_vector)
extroversion_alphas <- as.data.frame(extroversion_alphas)
extroversion_alphas <- rename(extroversion_alphas, c("country_vector"="country","ext_alpha_vector"="EXT_ALPHA"))

```

Table 1: Internal consistency of subjective well-being scale by country.

```{r}

country_vector <- c()
swl_alpha_vector <- c()

for (temp.country in levels(TEMP1$country)) {
  swlalpha <- subset(TEMP1, country==temp.country, swl1:swl5)
  a<-alpha(swlalpha, keys = c(1,1,1,1,1))
  alpha_output<-as.numeric(a$total[1])
  country_vector<-append(country_vector, temp.country)
  swl_alpha_vector<-append(swl_alpha_vector,alpha_output)
  }

wellbeing_alphas <- cbind(country_vector, swl_alpha_vector)
wellbeing_alphas <- as.data.frame(wellbeing_alphas)
wellbeing_alphas <- rename(wellbeing_alphas, c("country_vector"="country","swl_alpha_vector"="SWL_ALPHA"))

```

Table 1: Correlations, CI's and significance tests

```{r}

#Exroversion ~ Satisfaction with Life
swl <- ddply(TEMP1, "country", summarise, corr=cor.test(ext_tot, swl_tot))
swl_df <- swl[seq(2, nrow(swl), by=9), ] 
swl_cor <- swl[seq(4, nrow(swl), by=9), ] 
swl_ci <- swl[seq(9, nrow(swl), by=9), ] 
swl_p <- swl[seq(3, nrow(swl), by=9), ] 

swl1 <- merge(swl_df, swl_cor, by="country")
swl2 <- merge(swl1, swl_ci, by="country")
swl2 <- rename(swl2, c("corr.x"="df","corr.y"="swl_cor","corr"="swl_ci"))
swl_tot <- merge(swl2, swl_p, by="country") 
swl_tot <- rename(swl_tot, c("corr"="swl_p-value"))

```

```{r}

extroversion_alphas
wellbeing_alphas


swl_tot
country.means_sds

#Merge datasets 
means_and_corrs <- merge(country.means_sds, swl_tot, by="country")
all_alphas <- merge(extroversion_alphas, wellbeing_alphas, by="country")
table1final <- merge(means_and_corrs, all_alphas, by="country")

#Convert 'df' variable to 'n', and then delete 'df' 
table1final$df <- as.numeric(table1final$df)
table1final$N <- table1final$df +2
table1final$df <- NULL

#Rename columns 
table1final <- rename(table1final, c("swl_p-value"="SWL_p"))

#Reorder items in table. 
table1final <- table1final[c("country", "N", "EXT_M", "EXT_SD", "EXT_ALPHA", "SWL_M","SWL_SD", "SWL_ALPHA", "swl_cor", "swl_ci", "SWL_p")]

names(table1final)

#Specify data type in each column 
table1final$country <- as.character(table1final$country)
table1final$N <- as.character(table1final$N)
table1final$EXT_ALPHA <- as.character(table1final$EXT_ALPHA) 
table1final$SWL_ALPHA <- as.character(table1final$SWL_ALPHA) 
table1final$swl_cor <- as.numeric(table1final$swl_cor)
table1final$swl_ci <- as.character(table1final$swl_ci)
table1final$SWL_p <- as.numeric(table1final$SWL_p)

#Export Table
write.csv(table1final, "Table 1.csv")

``` 

NoTE: Final table tidied in Microsoft Excel 2013, and then transferred to manuscript. 